#include "PO_MoM.h"
#include "SolveEquation.h"
using namespace blitz;

void GetI()
{
	I.resize(iWire);
	double eps = 0.001;
	CG(iWire, V, 1000, eps);
	for(int i=0;i<iWire; i++)
	{
		std::cout<<abs(I(i))<<std::endl;
	}
}
void CG(int iUnknown, Array<dComplex,1> &zCb, int iItmax, double eps)
{
	Array<dComplex,1> zCtemp(iUnknown), zCx(iUnknown), zCr(iUnknown), zCd(iUnknown);

	zCx = 0;

	firstIndex i;
	secondIndex j;

	zCr = zCb - sum(Z(i,j) * zCx(j),j);

	zCd = sum(conj(Z(j,i)) * zCr(j),j);

	zCd = zCd / sum(abs(zCd) * abs(zCd));

	int k=1;
	while (k++<iItmax)
	{
		zCtemp = sum(Z(i,j) * zCd(j),j);
		zCx += zCd / sum(abs(zCtemp) * abs(zCtemp));
		zCr -= zCtemp / sum(abs(zCtemp) * abs(zCtemp));
		zCtemp = sum(conj(Z(j,i)) * zCr(j),j);
		zCd += zCtemp / sum(abs(zCtemp) *abs(zCtemp));
		double Eps0 =  sqrt(sum(abs(zCr) * abs(zCr))) / sum(abs(zCb) *abs(zCb));
		std::cout<<k-1<<"		"<<Eps0<<std::endl;
		if (Eps0<eps) break;
	}
	I = zCx;
}